
! vahovaci funkce pro robustni odhady statistickych parametru

module cutfun

  real,parameter,private :: huber_a = 1.345
  real,parameter,private :: hampel_a = 1.7, hampel_b = 3.4, hampel_c = 8.5
  real,parameter,private :: Andrews_a = 2.1, Pi = 3.1415929
  real,parameter,private :: Tukey_a = 6.0

contains 

! -------------------------------------------------------------
! 
! Huberovova funkce

function huber(x)

  real :: x,huber
 
  if( abs(x) < huber_a )then
     huber = x
  elseif( abs(x) >= huber_a )then
     huber = sign(huber_a,x)
  else
     huber = sign(huber_a,x)
!     write(*,*) x,huber_a   ! x = NaN
!     stop '1'
  endif

end function huber

function dhuber(x)

  real :: x,dhuber
  
  if( abs(x) < huber_a )then ! absolutni hodnota z derivace x
      dhuber = 1.0
  elseif( abs(x) >= huber_a )then
     dhuber = 0.0
  else
     dhuber = 0.0
!     write(*,*) x,huber_a
!     stop '2'
  endif

end function dhuber
 
!-----------------------------------------------------------------------
!
! Hampelova funkce

function hampel(x)

  real :: x,hampel
 
  if( abs(x) < hampel_a )then
     hampel = x
  elseif( hampel_a <= abs(x) .and. abs(x) < hampel_b )then
     hampel = sign(hampel_a,x)
  elseif( hampel_b <= abs(x) .and. abs(x) < hampel_c )then
     hampel = hampel_a*(x - sign(hampel_c,x))/(hampel_b - hampel_c)
  elseif(  abs(x) >= hampel_c )then
     hampel = 0.0
  else
     hampel = 0.0
!     write(*,*) x,hampel_a
!     stop '3'
  endif

end function hampel

function dhampel(x)

  real :: x,dhampel
  
  if( abs(x) < hampel_a )then ! absolutni hodnota z derivace x
      dhampel = 1.0
  elseif( hampel_b <= abs(x) .and. abs(x) < hampel_c )then
     dhampel = - sign(hampel_a/(hampel_b - hampel_c),x)
  else
     dhampel = 0.0
  endif

end function dhampel

!-----------------------------------------------------------------------
!
! Andrewsova funkce

function Andrews(x)

  real :: x,Andrews
 
  if( abs(x) < Andrews_a*Pi )then
     Andrews = sin(x/Andrews_a)
  else
     Andrews = 0.0
  endif

end function Andrews

function dAndrews(x)

  real :: x,dAndrews
  
  if( abs(x) < Andrews_a*Pi )then ! absolutni hodnota z derivace x
      dAndrews = cos(x/Andrews_a)/Andrews_a
  else
     dAndrews = 0.0
  endif

end function dAndrews

!-----------------------------------------------------------------------
!
! Tukeyova funkce

function Tukey(x)

  real :: x,Tukey
 
  if( abs(x) < Tukey_a )then
     Tukey = x*(1.0 - (x/Tukey_a)**2)**2
  else
     Tukey = 0.0
  endif

end function Tukey

function dTukey(x)

  real :: x,dTukey
  
  if( abs(x) < Tukey_a )then ! absolutni hodnota z derivace x
      dTukey = (1.0 - (x/Tukey_a)**2)*(1.0 - 5.0*(x/Tukey_a)**2)
  else
     dTukey = 0.0
  endif

end function dTukey

end module


